#include<iostream>
#include<arm_neon.h>
#include<sys/time.h>
#include <ratio>
#include<pthread.h>
#include<semaphore.h>
#include<semaphore.h>
using namespace std;
#define N 4800
#define num_thread 8
int n=100;
float m[N][N],init_m[N][N];
void init(int n){
    for(int i=0;i<n;i++){
        for(int j=0;j<n;j++)
            init_m[i][j]=0;
        init_m[i][i]=1.0;
        for(int j=i+1;j<n;j++)
            init_m[i][j]=rand()%1000+1;
    }
    for(int k=0;k<n;k++)
        for(int i=k+1;i<n;i++)
            for(int j=0;j<n;j++)
                init_m[i][j]=int((init_m[i][j]+init_m[k][j]))%1000+1.0;
}
void reset(int n){
    for(int i=0;i<n;i++)
        for(int j=0;j<n;j++)
            m[i][j]=init_m[i][j];
}
void serial(int n){ //串行
    struct timeval h, t;
    float time = 0.0;
    gettimeofday(&h, NULL);
      for(int k=0;k<n;k++){
          for(int j=k+1;j<n;j++)
              m[k][j]/=m[k][k];
          m[k][k]=1.0;
          for(int i=k+1;i<n;i++){
              for(int j=k+1;j<n;j++)
                  m[i][j]-=m[i][k]*m[k][j];
              m[i][k]=0;
          }
      }
    gettimeofday(&t, NULL);
    time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"serial: "<<time<<" ms"<<endl;
    return;
}

typedef struct{
    int t_id;   //线程id
}threadparam_t;
sem_t sem_div;  //信号量定义
sem_t sem_elimination;
pthread_barrier_t barrier;

void *threadfunc(void *param){
    threadparam_t *p=(threadparam_t *)param;
    int t_id=p->t_id;
    for(int k=0;k<n;k++){   //除法操作
        if(t_id==0){
            float32x4_t vt=vld1q_dup_f32(m[k]+k);
           int j;
           for(j=k+1;j+4<=n;j+=4){
               float32x4_t va =vld1q_f32(m[k]+j);
               va=vdivq_f32(va,vt);
               vst1q_f32(m[k]+j,va);
           }
           for(;j<n;j++)
               m[k][j]=m[k][j]/m[k][k];
           m[k][k]=1.0; 
            }
        else sem_wait(&sem_div);    //阻塞，等待完成除法操作
        if(t_id==0){
            for(int i=1;i<num_thread;i++)
                sem_post(&sem_div);     //唤醒，准备进行消去操作
        }
        for(int i=k+1+t_id;i<n;i+=num_thread){  //处理消去操作时，每个线程以8行为间隔处理某一行
            float32x4_t vaik=vld1q_dup_f32(m[i]+k);
               int j;
               for(j=k+1;j+4<=n;j+=4){
                   float32x4_t vakj=vld1q_f32(m[k]+j);
                   float32x4_t vaij=vld1q_f32(m[i]+j);
                   float32x4_t vx= vmulq_f32(vaik,vakj);
                   vaij = vsubq_f32(vaij,vx);
                   vst1q_f32(m[i]+j,vaij); 
               }
               for(;j<n;j++)
                   m[i][j]=m[i][j]-m[i][k]*m[k][j];
               m[i][k]=0;
        }
        pthread_barrier_wait(&barrier);
    }   
    pthread_exit(NULL);
}

void pthread(){
    struct timeval h, t;
    float time = 0.0;
    gettimeofday(&h, NULL);

    sem_init(&sem_div,0,0);     //初始化信号量
    sem_init(&sem_elimination,0,0);
    pthread_barrier_init(&barrier,NULL,num_thread);
    //创建线程
    pthread_t handles[num_thread];
    threadparam_t param[num_thread];
    for(int t_id=0;t_id<num_thread;t_id++){
        param[t_id].t_id=t_id;
        pthread_create(&handles[t_id],NULL,threadfunc,(void*)&param[t_id]);
    }
    for(int i=0;i<num_thread;i++)
        pthread_join(handles[i],NULL);
    sem_destroy(&sem_div);
    sem_destroy(&sem_elimination);
    pthread_barrier_destroy(&barrier);
    
    gettimeofday(&t, NULL);
    time += 1000*(t.tv_sec - h.tv_sec) + 0.001*(t.tv_usec - h.tv_usec);
    cout<<"pthread: "<<time<<" ms"<<endl;
}

int main(){
    cout<<"n="<<n<<endl;
    init(n);
    reset(n);
    serial(n);
    reset(n);
    pthread();
    n=200;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pthread();
    n=400;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pthread();
    n=800;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pthread();
    n=1600;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pthread();
    n=3200;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pthread();
    n=4800;
    cout<<"n="<<n<<endl;
    reset(n);
    serial(n);
    reset(n);
    pthread();
    return 0;
}